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ABSTRACT 

The problem of solving perturbatively the equations describing the evolution of self- 
gravitating collisionless matter in an expanding universe considerably simplifies when 
directly formulated in terms of the gravitational and velocity potentials: the problem 
can be solved exactly, rather than approximately, even for cosmological models with 
arbitrary density parameter 51. The Eulerian approach we present here allows to calcu- 
late the higher-order moments of the initially Gaussian density and velocity fields: in 
particular, we compute the gravitationally induced skewness of the density and velocity- 
divergence fields for any value of f2, confirming the extremely weak J7-dcpendence of the 
skewness previously obtained via Lagrangian perturbation theory. Our results show that 
the separability assumption of higher-order Eulerian perturbative solutions is restricted 
to the Einstein-de Sitter case only. 
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1 INTRODUCTION 

Any model of large-scale structure formation based on the gravitational instability has to deal with the problem of the nonlinear 
evolution of the density and peculiar velocity fields. In the simplest of these scenarios, the mass content of the universe is 
treated as a cold fluid of dust. In spite of the simplicity of the dust fluid model, the analysis of the nonlinear gravitational 
evolution is extremely involved, the Poisson, continuity and Euler equations being highly nonlinear and nonlocal beyond the 
linear regime, even in their Newtonian version (Peebles 1980; see also Kofman & Pogosyan 1994, and references therein). 
In the attempt of addressing fundamental cosmological issues, a big effort has been devoted to the theoretical treatment of 
nonlinear gravity, both in the Eulerian and Lagrangian formulation (e.g. Shandarin & Zel'dovich 1989; Sahni & Coles 1994). 
The basic Eulerian approximation scheme is perturbation theory: the density contrast 8 = (p — pb)/pb (p being the density 
field and pi, its mean background value) and the peculiar velocity u are expanded about the background solution S = 0, u = 0, 
namely S = £ 5 (n) and u = £ u (n) with <5 (n) = 0(5 (1) ") and u (n) = 0(u (1)n ), where S w and u (1) correspond to the linear 
solutions, then the differential equations for any <r n ' and vS n ' are solved (Doroshkevich & Zel'dovich 1975; Peebles 1980; Fry 
1984; Goroff et al. 1986). Thus, for example, the question of whether power can be transferred from large to small scales 
has been examined in the framework of the second-order approximation by Juszkiewicz (1981), Vishniac (1983), Juszkiewicz, 
Sonoda & Barrow (1984), Coles (1990), Suto & Sasaki (1991), Makino, Sasaki & Suto (1992), Jain & Bertschinger (1994) and 
Baugh & Efstathiou (1994). The weakly nonlinear deformation of the original Gaussian density and velocity distributions in 
terms of the gravitationally induced higher-order moments - and the corresponding correlation hierarchy - is instead analysed 
by Peebles (1980), Fry (1984), Goroff et al. (1986), Grinstein et al. (1987), Grinstein & Wise (1987), Bernardeau (1992a,b), 
Juszkiewicz, Bouchet & Colombi (1993), Bernardeau (1994), Catelan & Moscardini (1994a,b), Juszkiewicz et al. (1994) and 
Lokas et al. (1994); extensions to more general intrinsically non-Gaussian models are discussed by Luo & Schramm (1993) 
and Fry & Scherrer (1994). Some of these results and related issues have been also successfully investigated in complementary 
N-body simulations by Hansel et al. (1985), Coles & Frenk (1991), Efstathiou et al. (1990), Efstathiou, Sutherland & Maddox 
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(1990), Moscardini et al. (1991), Messina et al. (1992), Weinberg & Cole (1992), Bouchet & Hernquist (1992), Lahav et 
al. (1993), Lucchin et al. (1994), Kofman et al. (1994) and Melott, Pellman & Shandarin (1994). 

Typically, with the exception of a few papers (Bernardeau 1992a, b; 1994), this kind of investigation has been confined 
within the limits of the Einstein-de Sitter cosmology, mainly because the condition fl = 1 enormously eases the solution of the 
dynamical equations for 5 and u: apart from theoretical expectations, there is however no conclusive observational evidence 
that our universe is actually flat (e.g. Peebles 1991; Coles & Ellis 1994). 

A powerful alternative approach is the Lagrangian description, where the individual fluid elements are followed during 
their motion. The most popular Lagrangian approach is the celebrated Zel'dovich approximation (Zel'dovich 1970a,b), widely 
used in cosmology, showing also extremely useful in reconstruction methods of initial conditions from velocity data (e.g. Nusser 
& Dekel 1992). 

A Lagrangian perturbative approach, pioneered by Buchert (1989), Moutarde et al. (1991) and Buchert (1992), is presently 
object of thorough investigation. The key difference with respect to the Eulerian approach is that one searches for solutions 
of perturbed trajectories about the linear (initial) particle displacement S 1 - 1 ' : S = ^ S*-™' where S'™' = 0(S ( - 1 -'), while the 
continuity equation for the density field is exactly satisfied. 

Solutions up to the third-order Lagrangian approximation have been obtained, in the case of an Einstein-de Sitter model 
by Buchert & Ehlers (1993) and Buchert (1994), and for more general Friedmann models by Bouchet et al. (1992) and Catelan 
(1995). Matarrese, Pantano & Saez (1994a,b) developed a relativistic Lagrangian treatment of the nonlinear dynamics of an 
irrotational collisionless fluid, which reduces to the standard Newtonian approach on sub-horizon scales, but is also suitable 
for the description of perturbations on super-horizon scales. 

The higher accuracy of Lagrangian perturbative methods, as compared to other currently studied approaches [such as 
the frozen flow (Matarrese et al. 1992) and frozen (or linear evolution of the) potential (Brainerd, Scherrer & Villumsen 1993; 
Bagla & Padmanabhan 1994) approximations, is discussed by Munshi & Starobinsky (1994), Bernardeau et al. (1994b) and 
Munshi, Sahni & Starobinsky (1994), still in the framework of an Einstein-de Sitter cosmology. Comparisons with N-body 
simulations in the fully developed nonlinear regime are given in Moutarde et al. (1991), Coles, Melott & Shandarin (1993), 
Melott et al. (1994) and Melott, Buchert & Weiss (1995). 

Finally, let us also mention some observational evidence that gravitational instability in the quasi-linear regime operates 
on large scales. Baumgart & Fry (1991) claim that the observed bispectrum obeys the hierarchical pattern; in Saunders et 
al. (1991), the distribution of QDOT-IRAS galaxies out to 140 ft, _1 Mpc (h being the Hubble constant in units of 100 km sec -1 
Mpc -1 ) is found positively skewed and inconsistent with the results of the standard cold dark matter model, even after a 
20 h~ Mpc smoothing (see also Coles & Frenk 1991; Park 1991); higher-order correlations and/or moments in the galaxy 
distribution have been recently also measured by Szapudi, Szalay & Boschan (1992), Meiksin, Szapudi & Szalay (1992), Bouchet 
et al. (1993), Gaztaiiaga (1992, 1994) and others, both in optically and IRAS selected catalogues; higher-order effects are 
detected on significantly large scales by the counts-in-cells analysis of N-body simulations of Baugh, Gaztaiiaga & Efstathiou 
(1994). The general observational evidence appears to be quite consistent with gravitational instability perturbation theory 
applied on random-phase initial conditions. 

The problem we deal with in this paper is thus of quite general interest. We present an Eulerian formalism which allows to 
solve in a rigorous perturbative way the equations governing a pressureless self-gravitating fluid in the expanding universe, with 
arbitrary density parameter Q. The problem notably simplifies when formulated directly in terms of the peculiar gravitational 
and velocity potentials. In other words, we look for perturbative solutions of the Bernoulli equation rather than the Eulcr 
equation: the Bernoulli equation is the corner-stone for several nonlinear approximation techniques (Nusser et al. 1991; 
Matarrese et al. 1992; Gramann 1993a; Lachieze-Rey 1993; Brainerd, Scherrer & Villumsen 1993; Bagla & Padmanabhan 
1994). 

More in detail, we derive a particular form of the equations governing the evolution of the gravitational potential ip and 
the velocity potential <E>: the dynamics of the fields ip and $ is determined by the longitudinal component of the current 
S\i carried by the fluctuation 8 with peculiar velocity u. An advantage of this treatment is that we can easily work out the 
higher-order velocity contributions u' n ' in a manifestly irrotational form in Eulerian configuration space [all the previous 
irrotational expressions for u' n ' were obtained in Fourier space (see Goroff et al. 1986)]. The importance of having explicit 
irrotational expressions for the velocity field at any perturbative order lies in the fact that a collisionless self-gravitating fluid 
cannot generate vorticity (before orbit mixing) and that conservation of circulation is applicable well beyond the linear regime 
(see the discussion in Bertschinger 1993). 

Apart from the evident mathematical simplification, there are several more reasons for working directly with the gravita- 
tional and velocity potentials. Unlike the density fluctuation field, the gravitational potential is required neither to have zero 
mean (since only its gradients are physically meaningful), nor to satisfy the positive-mass constraint, 5 > — 1; moreover, since 
it depends directly on the mass distribution, it is bias-independent. On the other hand, the velocity potential has been recently 
recognized as a powerful tool to test the cosmological power spectrum, becoming directly measurable from observational data 
(see Heavens 1993; Baugh & Efstathiou 1994). 
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There are two major previous attempts to solve the same problem we study in this work. Martel & Freudling (1991) 
applied Eulerian perturbation theory, but in the attempt to get second-order solutions for the density contrast and peculiar 
velocity fields in non-flat models, they forced the intrinsic non-separability of the Eulerian dynamical equations. In Bouchet et 
al. (1992), the separability is also assumed, but for the higher-order Lagrangian perturbative ansatz; the separability hypothesis 
is also implicit in Gramann (1993a,b), and Lachieze-Rey (1993). Actually, the general validity of the separability ansatz in 
Lagrangian perturbation theory has been recently proved by Ehlers & Buchert (1995, in preparation; see also Appendix B in 
Catelan 1995). 

As an application of our Eulerian formalism, we finally calculate the density skewness Szs induced by the gravitational 
evolution, starting from Gaussian initial conditions, as well as the skewness of the velocity-divergence, £3,,. 

Before concluding this Introduction we would like to mention an important difference between the Eulerian and Lagrangian 
perturbative approaches, which is relevant for the present work. Eulerian perturbation theory, at any order, cannot lead to 
shell-crossing and consequent vorticity generation by multi-streaming, contrary to the Lagrangian perturbation theory, and, 
what is most important, contrary to the correct evolution which does lead to caustic formation. This is because the perturbative 
Eulerian velocity field, in a given position x and at time t, is obtained at any order by a set of (generally nonlocal) quantities 
which are functions of the position x itself (and time): no memory of the initial (Lagrangian) position q of the particle which 
is passing there at that time can be present; in other terms, the system does not contain particle inertia. Clear evidence of this 
fact is provided by the frozen flow approximation (Matarrese et al. 1992), which precisely exploits this property of the linear 
Eulerian velocity field (but the same would be true at higher orders). Therefore, provided the initial velocity field is sufficiently 
smooth, which will be always the case after suitable low-pass filtering, no caustic ever forms, as this would correspond to a 
discontinuity in the velocity field. Of course, one might ask whether this is a real weakness of Eulerian perturbation theory 
itself, compared to the Lagrangian one. Certainly, it is a limitation if one is interested in using its results to reconstruct the 
exact position of each particle after nonlinear evolution (e.g. by extending the method of the frozen flow approximation to 
higher orders). However, one is often interested in knowing only the coarse-grained features of the density and velocity fields 
in physical space: whether Eulerian theory at some perturbative order can reproduce such features remains in our opinion an 
open question. For the purpose of the present work the absence of multi-streaming is a clear advantage, for it allows the use 
of the velocity potential at any time during the nonlinear evolution. 

The plan of the paper is as follows. In Section 2 the Hamilton-Jacobi equation and the fluid evolution equations are 
derived starting from the single particle Lagrangian in a generic Friedmann universe; the Hamilton-Jacobi equation allows an 
alternative derivation of various approximation schemes, such as the frozen flow and Zel'dovich algorithms, which is sketched 
in Appendix Al, and shows an important analogy with the Bernoulli equation for the velocity potential of the fluid. In Section 
3, the second-order solutions for the gravitational and velocity potentials are worked out. These results allow to calculate the 
gravitationally induced density and velocity-divergence skewness and their dependence on the density parameter: in Section 4 
we compare these results with the corresponding ones obtained so far by means of more complicated approximation techniques 
or applying Lagrangian perturbation theory and transforming back to real space. Conclusions are drawn in Section 5. Technical 
appendices are also given. 



2 DYNAMICS OF SELF-GRAVITATING COLLISIONLESS MATTER IN NON-FLAT UNIVERSES 

We assume that, at the era of large-scale structure formation, the mass content of the universe is in the form of a collisionless 
and irrotational fluid with negligible velocity dispersion. Let x be Eulerian comoving spatial coordinates; physical distances 
are obtained by r = a(t)x, where a(t) is the scale factor and t the cosmic time. For an Einstein-de Sitter universe (density 
parameter Q — 1 and vanishing cosmological constant), a(t) oc t 2 ^ A . However, in what follows, we shall consider more general 
non-flat Friedmann models. In the general case a suitable time coordinate is the variable D(t), in that the dynamical equations 
assume very simple forms: D(t) is the growing mode of linear density fluctuations, whose dependence on t is in general quite 
complicated (e.g. Peebles 1980); with whole generality D(0) = at the initial time; in the flat case, D(t) oc a(t) oc t 2 ' [i . 

In this section we rederive the dynamical equations for our pressureless self-gravitating fluid directly from the single 
particle Lagrangian in a non-flat model. The Hamilton-Jacobi formulation of the problem is shown to be strictly related to 
the Bernoulli equation for the velocity potential; in terms of it the Zel'dovich and frozen flow approximations may be simply 
reformulated, as described in Appendix Al. 



2.1 From the particles to the fluid description 

Let Ct be the Lagrangian for the motion of a particle with mass m (e.g. Peebles 1980) 

£i(x,x,() = -ma x — m0(x, t) , (1) 



4 P. Catelan, F. Lucchin, S. Matarrese and L. Moscardini 

where 4> is the gravitational potential, related through the Poisson equation to the mass fluctuations 8, \7 2 cj> — 4-kG pb(i 2 8. 
The action of the particle reads S = J dtCt —J dD Co, where the last equality is justified by the fact that S is a scalar. Here 
dD = D(t)dt; in what follows, a dot indicates the operator d/dt, and a prime d/dD. From now on, for simplicity, we shall 
write C = Cd = D~ 1 Ct, and £(x,x',D) = ^rna 2 Dx' 2 — mD^ 1 (f>(x, D). An alternative form of the Lagrangian C may be 
obtained by suitably rescaling the gravitational potential cf> and the particle mass m. Indeed, introducing the two quantities 



ip = [2D/3e(n)a 2 D 2 ] <f> and fj,(D) = ab/y/D/e(Q) 



m, after some algebra one gets the final expression 



c ^' D) -^wmw^- 3 ^ D) }- (2) 

The function e(fi) is defined in terms of the logarithmic rate of growth of mass fluctuations, /(fi) = dlnD/d\na, through 
e(fi) = n/[f(il)] 2 (Gramann 1993a,b). Note that e(fi) depends weakly on Q, since e(fi) « fT 1/5 , once /(fi) w fi 3/5 is used 
(see Peebles 1980). The Einstein-de Sitter case is easily recovered with the replacements D — > a, e — > 1 and /i — > |o t^" 1 m 
[here to is an arbitrary reference time and ao = o(io)]- In terms of the potential tp the Poisson equation reduces to 

D V V(x, -D) = 5(x, D) . (3) 

Instead of using the velocity x' as independent variable, one can alternatively introduce the conjugate momentum p = M£r = 
aWD/eJJf) fix' . In such a case one adopts the Hamiltonian point of view, where the Hamiltonian reads 

«(x,p,I>) = p.x'-/:(x,x' ) I>) = f +1 — ^ y(x,£>) . (4) 

2a t i v /D/e(n) 2 ^D/e(Q) 

By varying the action S with respect to x and x' one obtains the Euler-Lagrange equations 

x" + ^(x' + V^)=0. (5) 

The three second-order Euler-Lagrange equations above are equivalent to the set of six first-order Hamilton equations, namely 
, _ &H _ p , _ &H _ 3 a/i 



dp au^/We ' '"' dx 



=V<^ . (6) 



Alternatively, one can use the Hamilton- Jacobi formulation, according to which the momentum p is written in terms of 
an action functional S, as p = VS, the latter obeying the Hamilton-Jacobi equation S^ +H(x, V5) = 0. It can be, however, 
useful to define a scaled action functional $5 through <S(x, D) = a-J D/e(Q.) fi $s(x, D) so that the Hamilton-Jacobi equation 
takes the simple form, 

ft + f<v*rf + 2fW,) = o. (7, 

So far we have only dealt with the dynamics of single particles. Our problem is, however, the description of an infinite set 
of particles interacting only through gravity. The standard description of this system is based on the Vlasov equation (or colli- 
sionless Boltzmann equation), which governs the evolution of the (comoving) one-particle distribution function /(x, p, D)dxdp, 
namely the probability of finding one particle in the infinitesimal phase-space element dxdp. 

By the Liouville theorem the distribution function is conserved along particle trajectories in phase-space, thus leading to 
the Vlasov equation, which in our case reads 

21 + , p -y/-g , a » S7,..gJ.=n (8) 

dD any/D/e(Q) 2 ^D/e(n) dp 

This equation, together with the definition of comoving mass density 

£>(x, D) = m n(x, D) = po(l + 8(x, D)) = m j dp f(x, p, D) (9) 

(here n is the comoving particle number density and po the mass density at to) and the Poisson equation (3), completes 
the description. It is however useful to attempt a fluid description of the system, which is achieved by introducing the first 
moment of the distribution function, namely the local streaming velocity, i.e. the mean velocity of the patch of fluid around x 

< m 1 frf PP /(x,p,D) 

u(x, D) — , ^ — . (10) 

ati^/D/e(n) Jdpf(x,p,D) 

From the Vlasov equation, we then obtain the continuity equation 
Ig+V-^uHO (11) 

and the momentum conservation equation 
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-j— +U0dpu a + o a ip = — opyQirctp) 

(where d a = d/dx a ; sum over repeated indices is understood), having defined the velocity dispersion tensor, 
1 J dpp a pi3f(x,p,D) 



*0 



aV-D/e(n) fdpf(x,p,D) 



UaU/3 



(12) 



(13) 



For vanishing initial velocity dispersion and vorticity, and as long as the system remains in the laminar regime we can 
write the distribution function in the single-stream form 

/(x, p, D) = n(x, D) Sd (p - a^/D/e{Q.) n V x $s(x, D) 



(14) 

[see also (Vergassola et al. 1994) for a related expression], with Sd the Dirac delta function, which immediately leads to 
the vanishing of the velocity dispersion tensor and to the validity of the fluid picture. The resulting form of the momentum 
conservation equation is indeed the Euler equation for a pressureless fluid in a general non-fiat universe (e.g. Shandarin 1994), 



du , „. 3e(fi) . „ „ 

_ + ( u .V)u + ^_2(u + V^)=0 



(15) 



Given the vanishing of the initial vorticity, Kelvin circulation theorem implies that no vorticity will be generated at later 
times, except for regions of multi-streaming (e.g. Doroshkevich 1973). In these conditions we can define a velocity potential 
$ via u = V4>. Replacing this into the Euler equation one then finally gets the Bernoulli equation 

d<S> 1 
dD + 2 ( 



+ , ( v*f + 3 ^ 



2D 



(<3> + <p)=0 



(16) 



which is formally identical to the scaled Hamilton- Jacobi equation (7). One should, however, keep in mind that the action 
$5 and the velocity potential $ generally represent different physical quantities. Moreover, while the validity of equation (7) 
has no restrictions at all, the Bernoulli equation (16) is only valid provided that i) one has vanishing initial vorticity and ii) 
vanishing initial velocity dispersion, and iii) the system is in the single stream (or laminar) regime, i.e. no orbit crossing has 
yet taken place. 



2.2 Dynamics of the gravitational and velocity potentials 

By expressing the density fluctuation S and velocity field u in terms of the gravitational and velocity potential, respectively, 
we can recast the continuity equation in the form (see also Kofman 1991) 

^ (DVp) + (1 + DV V) V$] = , (17) 

whose general solution reads 

^(£>Vv>) + (l + .DVV)V$ = VAT. (18) 

A standard way to single out the irrotational part of the latter equation is the following: according to Helmholtz's theorem 
(e.g. Morse & Feshbach 1953), the last term in the l.h.s may be written as (up to an irrelevant harmonic function) 



(DV 2 <p)V$ = <5u = VF + V A T , 

in such a way that the solenoidal vector T in equation (18) is automatically canceled out: 

^(DV^)+V$ + VF = 0. 

The function F just introduced has to satisfy the equation 

V 2 F = V ■ (<5u) = <5V 2 $ + V<5 • V<& . 

To summarize, one can write the set of equations for the cosmological potentials as follows: 
9 



(19) 



(20) 



(21) 



aD] W*D) 



+ $(x,D)+F(x,D) =0 
3e(fi) 



V$(x,D) 



+ 



2D 



$(x, D) + <p(x, D) 



0, 



(22) 



V 2 F(x, D) = DV- [VV(x, D) V$(x, £>)] 
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Some remarks are appropriate. The auxiliary potential F - which closes our system of scalar equations - is related to the 
flux of matter originated by the mass current j = <5u carried by the density fluctuation S moving with peculiar velocity u: 
only the longitudinal component j|| determines the dynamics [a related quantity is introduced by Bertschinger & Hamilton 
(1994), although in a different context]. Furthermore, the function F is, by construction, at least a second-order quantity in 
the fields S and u, because j is so. This fact allows us to easily linearize the equations: the first-order form of the equations 
(22) is 

<,(i) + S<i)=0, 



dip 



(i) 



(23) 



dD 



= 



so that the linear regime is described in terms of a single potential, say tp (see Kofman 1991). Note that the second of these 

(x) and thus $(x, D) = <J> (1) (x) = 



(i) 



equations implies that the gravitational potential is time-independent, i.e. <fi(x, D) = <p 
— (^'(x), in the linear regime. 

Only a few exact solutions of the set of equations (22) are known, mostly endowed with some symmetries, i.e. with 
restrictions in the initial conditions. A possible alternative strategy is to seek perturbative solutions: this program will be 
carried out in the next section. 



3 EULERIAN PERTURBATION THEORY: SECOND-ORDER APPROXIMATION 



Let us now solve the previous equations for the potentials according to the following general Eulerian perturbative ansatz: 

V >(x,D) = l p {1) (x) + ^ 2) (x,D) + --- = lfl {x) + D^ 2 (x, £>) + ••■ , (24) 

$(x,D)=$ (1) (x) + $ (2) (x,D) + --- = $i(x) + D$ 2 (x, J D) + --- , (25) 

explicitly up to second-order: of course, the expansion may be continued to higher-order terms. The functions ip^'(x) and 
<lr '(x) are the linear potentials, as already discussed. Note that in our notation, for instance, (r (x, D) = D V ip^ (x), but 
<5i(x) = VVi(x), where S (1) (x,D) = D6i(x); instead: u (1) (x) = V$ (1) (x) = V$i(x) = ui(x). Similarly for the higher-order 
terms: the expansions (24) and (25) obviously generate analogous expansions for the corresponding fields S and u. 

We assume that the general ansatz in (24) and (25) actually catches all the physical information contained in perturbation 
theory: in this sense we stress that, in particular, the second-order solutions <^ 2 '(x, D) and $"'(x, D) are non-separable in 
the variables x and D. On the contrary, what is usually done is to assume that all the temporal information contained e.g. 
in <p( 2 '(x, D) can be factored out in the form ip( 2 '(x, D) — D <p 2 (x): this is not allowed a priori, it applies in the case of the 
Einstein-de Sitter universe, but not in a generic Friedmann model. To be convinced of the previous statement it is enough 
to suppose, ab absurdo, that separability actually holds, then show that such an ansatz is not a solution of the fundamental 
equations; we give this proof in Appendix A2. 

In terms of the perturbed quantities, our equations become 

d 



D jD' fi2 ^' D) + 2 ^(*' D ) + *2(x, D) + F 2 (x) 







Dq B ^{x,D) + 



1 + 



3e(n) 



$ 2 (x,D) + ^gW,D) + i 



V$i(x) 



0, 



(26) 



V 2 F 2 (x) =V- [vV(x)V$i(x)] = V' 



<5i(x)ui(x) 



where the dependence on the variables x and D is explicitly shown and, after noting that -F(x, D) is actually a second- 
order quantity, we have defined the scaled function F?,{x) = D^ 1 F^ 2 ' (x, D); the solution of the last equation above can be 
immediately obtained as a convolution in Fourier space, F 2 (k) = J rfxF 2 (x) c lk x : 



F 2 (k) 



/,;-' 



dkidk 2 

(2tt) 3 



fo(ki + k 2 -k) 



1 + 



kik 2 



"2 



fc 2 <^i(ki) fc|<^i(k 2 ) 



(27) 



where e.g. k — |k| : F2 is determined by the initial condition on tpi, while the kernel 1 + ki • k 2 /fc 2 describes the effects of the 
nonlinear evolution. 

The differential equations in (26) are coupled in the fields <p 2 (x, D) and $ 2 (x, D); to decouple them, we differentiate the 
first once, to get 



_ d v? 2 d^2 , , _ 



(28) 
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By substituting into the second equation of the set (26) we obtain 



D 



d 2 ip 2 
dD 2 



4 + 



3e(fi) 



D 



dtp 2 
dD 



3e(n) 



H>2 = -^{Vipif 



1 



3e(n) 



F 2 



(29) 



which turns out to be an equation for the single potential y> 2 . Now it should be noticed that, while </5 2 (x, D) is not separable 
in space and time, the sum 



V>(x, J D) = ^ 2 (x,D)+F 2 (x) 

is exactly separable; in fact, after substitution in equation (29), we find 



D 2^± 

dD 2 

where 



<7(x) = - 



3e(fi) 



Vv?i(x) 






F 2 (x) 



2 + 



3e(fi) 



V = C(x) 



Defining at this point the function 
A{0) = 2 + | e(n) = A{D) , 

we can summarize our results in terms of the following set of equations, 
d 2 / .,„A„ 9 



(30) 



(31) 



(32) 



(33) 



D' 



<9D 2 



2 + A(D))Z3— +A(£>) 



^(x,I>) = C(x) 



<p 2 (x,L>) =V(x,D)-F 2 (x) 



f) 



(34) 



$ 2 (x, D) = -D— V2 (x, D) - 2^ 2 (x, D) - F 2 (x) , 



which exactly defines the second-order potentials <£ 2 (x, D) and "l? 2 (x, D). In particular, since F 2 (x) is determined by the initial 
conditions, the last equation of the set (34) gives the velocity potential $ 2 (x, D) once <^ 2 (x, D) is known: the main task is to 
solve the first of the equations in (34); we address this problem in the next subsection. 



3.1 Second-order perturbative solutions 

We have previously shown that, although the second-order potential <y3 2 (x, D) is non-separable, the function <p 2 (x, D) + F 2 (x) 
is so: since the r.h.s. of the first equation in (34) does not depend on D, we can assume (up to an additive constant, which 
can be always set to zero) 

ip(x,D) = B(D)C(x), (35) 

where the function B(D) is a function of the only variable D. Therefore, the second-order potential <p2 (x, D) takes the form 



V? 2 (x,D) 



l-B(D) 



F a {x) + -B(D) 



Vpi(x) 



The second-order velocity potential then reads 



$ 2 (x,D) 



1 - 2B(D) - DB'(D) 



F 2 (x)-i 



2B{D)+DB'(D) 



Vpi(x) 



The last relation shows that the velocity u may be indeed written as the gradient of a scalar function in {x}-space: 



u(x,D) = V 



$ 



(i) 



1 - 2B(D) - DB'(D) JF (2) (x,D) - ^D[2B{D) + DB'(D))u <1)2 



(36) 



(37) 



(38) 



explicitly up to second-order; this expression of the velocity field is thus manifestly irrotational in physical space [see also 
Catelan & Moscardini (1994b); for its Fourier components see Goroff et al. (1986)]. 

We can now furthermore derive the final expressions of the second-order corrections <r 2 '(x, D) = D 2 5 2 (x, D) and V • 

u (2) (x,D) e])V-u 2 (x,D): 



<5 (2) (x,D) = 



1-5(13) 



<5 (1) (x,D) 2 - D (u (1) • V) <5 (1) (x,D) + D 2 B(D) ^ (cW 1 



)P\ 



(39) 



a/3 
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and 

- D V • u (2) (x, D) 



1- 25(D) - DB'(D) 5 (1) (x,D) 2 -D(u (1) • V)<5 (1) (x,D) 



+ D< 



2B(D) + DB'(D)] Y^ (9 Q u (1 



l)/3\ 



(40) 



ad 



Combining these solutions, we can also derive a relation between the density contrast 8 and the divergence of the velocity 
field V • u, namely 



DV -\i = 8~2D z 



DB{D) 



Mm) , 



(41) 



where /i2(ui) is the second-order invariant of the initial deformation tensor 9 a 9^$i(x), which is in general different from zero 
in three-dimensional systems: 2/42 (ui) = (V • ui) 2 — ^ g (<9 a wJ (x)J = 2 (A1A2 + A1A3 + A2A3). The quantities Ah, h = 1, 2, 3, 
are the eingenvalues of d a uf(x). The relation (41), which quantifies the second-order deviation of V • u from 5, may be 
compared with similar relations in the literature (e.g. Giavalisco et al. 1993; Mancinelli et al. 1993; Mancinelli & Yahil 1994). 
The importance of the shear term ^ . [d a ul (x)j in the density- velocity-divergence relation, even in the case of laminar 
flow, is also discussed by Mancinelli & Yahil (1994). 

The corresponding expressions of 8^ 2 ' and vr 2 ' for the flat limit case (see Peebles 1980, where a different gravitational 
potential is used; Goroff et al. 1986, where a different normalization of the velocity field is chosen; also Appendix A3 for an 
independent derivation) read: 



5 (2) (x, a) = I 8 W (x, a) 2 -a (u (1) • V) 8 {1) (x, a) + ^a 2 Y / (d a u w ?) 



and 



a V • u< 2 > (x, a) = 5 S V> (x, a) 2 - a (u« ■ V) <5 (1) (x, a ) + ^a 2 J2 (d a u^ j ) 



(42) 



(43) 



a a 



where, in the Einstein-de Sitter model, S^ 2 ' scales like a 2 and u' 2 ' like a, in that they are separable in space and time, 
<r 2 ' (x, a) — a 2 82 (x) and \v 2 ' (x, a) = a 112 (x) . 

The function B(D), defined in (35), allows to specify the (second-order) growth of the potentials <^(x, D) and $(x, D): it 
is a solution of the differential equation 

1 r ~ ■~ 1 ~ d - K ■ (44) 



LY 



dD 2 



2 + A(D) 



D—+A(D)}B(D) = 1 



The initial condition (at the beginning: D — 0) for the function B(D) may be easily found comparing equations (39) and 
(42): B(D = 0) = 2/7. Comparing instead equations (40) and (43), we get a relation for the derivative of B(D) at the initial 
time, namely DB'(D)\d=o — 0. 

Since the dependence of A(D) on D is definitively non-trivial, we transform to a new variable r which, e.g. in the case 
Q < 1, is defined by (see e.g. Peebles 1980) 

1 



n- 1 



(45) 



where 0<r<oo:r = corresponds to the Big Bang (where one asymptotically approaches the Einstein-de Sitter model). 
The previous differential equation can now be written as 



J(t 



dr 2 



J'(t) 



Mr] 



J(r)-^ + A(r)} B[r)= I 



(46) 



and the initial conditions keep the same form. Here and in what follows, a prime denotes differentiation with respect to r. To 
maintain the notation concise we have defined 



J(t) 



D(r) 
D'(t) 



^rfiry 1 



and J(t) — > t in the limit r — » 0. The growing mode D(t) reads (e.g. Peebles 1980): 



D(t) = 1 + - + 3/2 L(t) 



with L(t) = lnl y/1 + t — y/r); D(t) — * 2r/5 in the limit r — ► 0. The function A(r) is now simply given by 
3 



A(r) = 2 + 



(1 + t)/(t)' 



(47) 



(48) 



(49) 
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where finally 


/(r) = -| 


3 v /r(r + l) + (3 + 2r)L(r) 


(3 + r) ^Jt(t + 1) + (3 + 3r)L(r) 



(50) 



and /(0) = 1. 

The analytical solution of Eq.(46) for the open case, with growing mode initial conditions, B(t)\ t =o — 111 and tB'(t)\ t= o — 
0, can be guessed by noting that there exists a close relation between our Eq.(46) and Eq.(7) in Bouchet et al. (1992), such 



that B = i 



1 + 



gill 



where E(t) is also reported in eq.(70) below. Indeed, the function 

n 2\ 



B v = \-m 



9 r /l + T r , . i 

4^(7F 1 + V~^ L(r)+ 2 



1+1 , L(r) 



(51) 



solves the differential equation (46) with the appropriate initial conditions. This result also shows the equivalence between the 
Eulerian perturbation theory, which is in general characterized by non-separable perturbative solutions, and the Lagrangian 
description, whose higher-order solutions are intrinsically separable (see the discussion in Appendix B of Catelan 1995 for the 
second-order case; the separability of higher-order Lagrangian modes has been recently demonstrated by Ehlers & Buchert 
1995, in preparation; Buchert, private communication). 

Incidentally we give here also the corresponding expressions directly in terms of the variable f2, which can alternatively 
be chosen as time variable: in the order, 



J(Cl)-- 
D{Cl) 



-fi(l-fi)/(fi) 

1 



(1 - Q) 3 / 2 



[(l + 2fi)v / I _ n + 3QL(Q.)] , 



/(0) = -f 



3VT^TT + (2 + SI) L(fi) 



(1 + 2fi) VT^Tl + 3SL L(Sl) 



(52) 

(53) 

(54) 



where L(SY) = In 



[l-VT=Sl)/Vsl 



and L(l) = 0. There are two approximations of the function f(S~l) usually found in 

the literature: f{Sl) « i7 3/5 (Peebles 1980) and /(fi) ^ Q 4/7 (Fry 1985). Also, we give all the formulae for the closed model 
(fi > 1) in Appendix A4. 



4 APPLICATIONS: THE GRAVITATIONAL SKEWNESS 

The previous results may be used to compute the gravitationally induced skewness of the density fluctuation, S, and divergence 
of the velocity field, r\ = — D V • u. We assume that the primordial potential tp is Gaussian distributed; then, for the density 
contrast the first non-zero contribution to the skewness is 



<<5 3 }=3(<5 (1) V 2) > + 0(«5 (1)6 ) 



(55) 



It is well known that this kind of calculations is better performed in Fourier space, namely we need to compute the term 
(8 w (k 1 ,D)8 w (k 2 ,D)8 (2) (k 3 ,D)}. From equation (39), one finally gets 



&(k,D)= / *£**fo( kl 



f :1 : , k 2 -k)J s (2) (ki,k 2 ;D)5i(ki)5i(k 2 

where the symmetrized kernel Ji is defined by the relation 



Jj 2) (ki,k 2 ;£>) 



l-B(D) 



1 /fci fo 

2 U 2 fei 



kl ' k2 +B(D)^ kl ' k2X2 



fcl fc2 V fcl ^2 

The gravitational skewness of the unfiltered density field may thus be written in the form 
dkdk' 



{51 = 6 



(2yr)t 



4 2) (k,k';D)P(k)P(k') 



(56) 



(57) 



(58) 



where P(k) is the primordial power spectrum of the density field. Performing the integration, we obtain the following expression 
for the density skewness parameter S3S: 

(5 A ) _ 3(5 (1)2 ,5 (2) ) _ 



S^{t) = 



(52)2 



(<5(!) 2 ) 2 



4B(r) , 



(59) 
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having emphasized the dependence on, e.g., the variable r. Note that, at the beginning B(0) — 2/7, and therefore Ss(0) — 34/7, 
which is the well-known value of the unfiltered density skewness in the Einstein-de Sitter model (Peebles 1980). 
Similar calculations can be done for the velocity-divergence field r\. It turns out that 

dki<ik 2 



m(KD) 



(2tt) 3 



S D (ki + k 2 - k) Ki 2) (ki , k 2 ; D) Si (ki ) Si (k a 



where the kernel 7-Q is defined by 



(2) 



1 - 2B(D) - DB'(D) 



tfj a) (ki,k 2 ;2J) = 

Then 



<^ 3 )=6 / ^-Ki 2 Hk,\,';D)P(k)P(k') 



1 /fci fa\ ki ■ k 2 



2 \ fc 2 fci / fci fc 2 



2B(D) + DB'(D) 



fci fc 2 



(2tt)« 



from which 



-^37, (r) = 



(v 3 ) 

iv 2 ) 2 



= 6 - 8B(t) - 4J(t)B'(t) 



(60) 



(61) 



(62) 



(63) 



Note that, for the Einstein-de Sitter model, we recover the standard result Ss v (0) = 26/7 (see Bernardeau 1994). 

4.1 Previous determinations of the skewness 

Both the results (59) and (63) can be directly compared with the corresponding results in the literature. 

4.1.1 Martel & Freudling (1991) 

In Martel & Freudling (1991), the second-order Eulerian equation for the density field is reduced to a separable version of it, by 
replacing Q 1 ' 2 — > ft in one of the terms. Martel and Freudling call this "almost-separability" of the perturbation equation for 
the density field, and the consequent approximation "MF approximation" (see also Martel 1995). In practice, this corresponds 
to assume that the expansion of the density contrast is of the factorized form 



<J(x,t) =D(t)6i(x) 



1 + k(t)C(x) + • 



up to the second-order term. The function C(x) is given by the expression 

1 ,„ „ , m2 



C(x) = §* 1 (x)- s ^V* 1 (x).V X (x) + 



167T 2 (5l(x] 



(d a d 0X (x)Y 



where X = J ^ x ' ^i( x ')/l x ' — x li furthermore, the function k(t) is a solution of the differential equation 



rV(T) + 2(n 3/5 + l 



^jtk' = (r + l)- 1 D(r) 



(64) 



(65) 



(66) 



with k(0) = k'(0) = and k(t) — * 2D(r)/7 in the limit r — » 0. It turns out that the density skewness is given by the 
expression 

K(T) 



S 3S (t) = 17 



D(r) 



(67) 



In the mentioned references no explicit estimate of Ss v is reported. 



4.1.2 Bouchet et al. (1992) 

The dynamical equations for the displacement S(q, r) from the Lagrangian (initial) positions q of the fluid elements are solved 
in Bouchet et al. (1992) according to the perturbative expansion 

S(q, r) = D(t) S (1) (q) + E(r) S (2) (q) + • • • . (68) 

up to the second-order correction [see also Catelan (1995) for the notation]. After (perturbativcly) inverting the relation 
x = q + 5(q, r), to get the Eulerian location of the mass elements, one obtains the following expression for the second-order 
density contrast <r 2 ': 



5^(x,r) = 



E/D 2 



D 2 



8 (1) (x, r) 2 - D (u (1) ■ V) S m (x, r) + ±L 1 + E / D 2 ^ (d a u {1 



l)/3\ 



(69) 



a/3 
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where the function E(j) is given by the relation 
EM = 




1 + T , L(t) 



(70) 



which, in the limit r — * 0, yields E — + —3D /7. The formal analogy between Eq.(69) and our expression Eq.(39), suggests 



that the two results coincide provided B 
results 

S m (t)=4-2 SW 



1 + 



Mil 
7577F 



which is indeed the case, as noticed above. The density skewness 

(71) 



D{rY ' 
while the skewness of the velocity-divergence reads (Bernardeau et al. 1994a): 

The divergence of the velocity field 6 in Bernardeau et al. (1994a) is defined in such a way that 6 = —f(Q)r] and thus 

s 3e = -/(a)- 1 ^. 



4.1.3 Bernardeau (1994) 

The Eulerian perturbative expansion here is based on the spherical model; it reads 

<S(x,r) = D(t)&(x) + i J D 2 (r)«5 1 (x) 2 + ■ ■ ■ , 

explicitly up to the second-order contributions. The function Di{t) is a solution of the differential equation 

3 3 



-D 2 



2 t(3 + 4t) , 

2 2(1 + t) 

and, in the limit r — > 

and the velocity-divergence one is 



-D z 



8 2 r-,/2 

3 rC 



2(1 + t) * 1 + T 
0, D 2 -> 34D 2 /21. The density skewness is 



S, 



> + 3- 



D' 2 (t) 



(73) 



(74) 



(75) 



(76) 



D(t)D'(t) ' 

Again, the divergence of the velocity field 9 in Bernardeau (1994) is such that 6 — —f(Q)r/. These Eulerian results for the 
skewness are identical to those obtained with Lagrangian methods by Bouchet et al. (1992): the expressions of S35 and S3,,, re- 
spectively in equations (71) and (72) and equations (75) and (76), coincide, in that the relation D 2 = |(2D 2 — E) actually holds. 

From this brief review we understand that, including ours, there are in the literature four independent calculations of 
Sz&{£l) and three of Sz n (£i). In particular, comparing equation (59) with equations (67) and (71), the two different expressions 



of B(t) have been used: i) B\(r) 



17 



«(t) 

D(t) 



is the equivalent of our function B(t) in Martel & Freudling (1991); ii) 



B2{T) = \ 



1 + 



E(t) 
OFF 



, in Bouchet et al. (1992), coinciding with our B(t), solution of Eq.(46). The disagreement with the 
Martel & Freudling results is summarized in Fig. 1. 

It should however be noticed that, in the observationally relevant Q range, the two determinations of the skewness differ 
only by a few percents, while the deviations become more conspicuous for low fl values. 



5 SUMMARY AND CONCLUSIONS 

In this paper we developed an Eulerian perturbative formalism which allows to solve the Newtonian equations governing the 
dynamics of a cold irrotational fluid in a generic Friedmann universe, with vanishing cosmological constant. In particular, the 
problem is notably simplified once described in terms of the gravitational and velocity potentials, respectively tp and <&. The 
general set of scalar equations for these potentials was obtained and it was shown that, in order to close this set one needs to 
introduce an auxiliary scalar potential related to the current j = <5u, carried by the density contrast S, moving with peculiar 
velocity u. 

Once the fundamental equations (22), which are nonlinear and nonlocal in the potentials, were obtained, the next step 
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Figure 1. f2-dcpcndcnce of the density skewness S35. The solid line represents our solution of eq.(59), which is equivalent to that of 
Bouchet et al. (1992) and Bcrnardcau (1994). The dotted line shows the corresponding expression obtained by Martcl & Frcudling (1991). 

was to solve them for arbitrary density parameter fi. We solved them according to the general Eulerian perturbative ansatz 
given in (24) and (25): apart from the approximation intrinsic in any perturbative approach, no further approximation 
was introduced. In particular, we focused in seeking the second-order perturbative solutions ip^ and $^ 2 ', whose explicit 
expressions are reported in equations (36) and (37), which we demonstrated not to be separable with respect to their space 
and time dependence; the separability assumption is correct only in the case of the Einstein-de Sitter universe, but not for 
more general Friedmann models. 

As a first application of our results, we have computed the second-order deviation of the velocity-divergence, V • u, from 
the density contrast, 8. We finally calculated the skewness parameters S35 and S3 V of the density fluctuation and velocity- 
divergence fields, respectively: these third-order moments quantify the asymmetric deformation of the original Gaussian density 
and velocity-divergence distributions, induced by gravity during the (weakly) nonlinear evolution. The weak dependence of 
S3S and S^,, on the density parameter Q, first demonstrated applying Lagrangian perturbation theory, was confirmed. 

Our results also clarify that the Eulerian perturbation theory, which is characterized by non-separable perturbative 
higher-order solutions (in the case of a general non flat Friedmann model), is equivalent to the corresponding Lagrangian 
description, whose perturbative solutions of the dynamical fluid equations are intrinsically separable with respect to the space 
and time coordinates. 

The present work provides a first attempt of building up a complete higher-order Eulerian perturbation theory: most 
previous papers only dealt with the construction of the density fluctuation and velocity-divergence fields [a remarkable 
exception being the work by Munshi & Starobinsky (1994), which however, besides adopting a different formalism is restricted 
to the flat case] . This kind of treatment should allow a direct comparison of the Eulerian approach with the Lagrangian one, 
which might lead to understanding many non-trivial aspects of the nonlinear dynamics of self-gravitating collisionless matter 
in the expanding universe. 
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APPENDIX Al: APPROXIMATE SOLUTIONS OF THE HAMILTON-JACOBI EQUATION 

A standard method to deal with the solution of the Hamilton- Jacobi equation (7) is through the simple ansatz 

$s = W(x) - £D , (Al) 

where W is Hamilton's principal function, and the constant £ plays the role of particle 'energy' per unit mass. In our case, 
however, this ansatz cannot work, since the peculiar gravitational potential, self-consistently determined from the Poisson 
equation (3), generally contains a non-trivial time dependence. Nevertheless, in most approximation schemes aimed at moving 
particles, the physical peculiar gravitational potential is practically replaced by some external, or 'mock', potential, e.g. 
derived from perturbation theory, which allows to simply close the fluid dynamical equations, without fully solving for the 
self-gravitation of the moving particles. Examples of this are: i) the Zel'dovich approximation (Zel'dovich 1970a, b), where one 
makes the ansatz <p(x,D) = — $(x, D), and ii) the frozen flow approximation (Matarrese et al. 1992), where one can write 
<p(x,D) = ipt 1 ' (x) — 3~jm \^ l P ( x )) [extending to Q 7^ 1 the expression in (Matarrese et al. 1992)]. Second-order type 
solutions can also be put in a similar form, but we will not deal with this issue here. 
Let us now analyze the two cases above in more detail. 



Al.l Zel'dovich approximation 

The action in this case has to satisfy the standard Hamilton- Jacobi equation for inertial motion, 

ff + i(V<M 2 =0, (A2) 

which can be immediately solved in the form (Al) with W(x) = — ^'(q) — V q y ( ' 1 '(q)(x — q) and £ — ^(VqV2 (q)) 2 . One 
then obtains the standard formula (e.g. Kofman 1991) 

* 5 (x,.D) = -v> (1) (q) + £ E ^, (A3) 

with x(q, D) = q - DV q <£ (1) (q). 



A1.2 Frozen flow approximation 

In this case, after replacing the appropriate expression for ip in the Hamilton- Jacobi equation, we immediately obtain 

$ s (x, J D) = - V 9 (1) (x), (A4) 

corresponding to vanishing particle 'energy', £ — 0. This conclusion is consistent with a simple picture of this approximation 
according to which the fluid elements behave as test particles, which just trace the initial stream-lines. 

The advantage of formulating these two approximations in terms of the Hamilton- Jacobi approach is that, once the 
solution is known, both mass and momentum conservation are exactly satisfied at all stages of the nonlinear evolution. 
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APPENDIX A2: FAILURE OF THE SEPARABILITY ANSATZ IN NON-FLAT MODELS 

In this appendix we explicitly demonstrate that, in the case of a generic Friedmann universe, any separable perturbative 
approximation is not a solution of the dynamical equations for the cosmological potentials tp and $. To avoid an exceeding 
proliferation of mathematical symbols, we use in this appendix the same notation of the main text. 

Let us therefore suppose that a possible class of solutions of the fundamental equations for the cosmological potentials 
may be written in the form 



S 2 ) 



<p(x, D) = <p w (x) + ip w (x, £>) + .•• = Vl (x) + D g v (D) ^(x) + • 



$(x, D) = $ w (x) + $ w (x, D) + 



(2), 



$i(x) + Dcm,(D)$ 2 (x) 



(A5) 
(A6) 



where the functions if2 (x) and $ 2 (x) are assumed to be independent of D, while the functions g<p{D) and g<s>(D) are assumed 
to depend only on D. The calculations performed in Section 2 can then be repeated up to Eqs.(36) and (37), which we rewrite 
in the form 

t2 



F 2 (x) + -B( J D) 



ff „(i% a (x) = - I -B(D) 

and 

<?i>(D)<& 2 (x) = l-2B{D)-DB'(D) 



V<pi(x) 



*b(x) - ^ 



2B(D)+DB'(D) Vv?i(x) 



(A7) 



(A8) 



where the two space-dependent functions <p2(x) and $ 2 (x) can be easily determined from the fi — > 1 limit (see Eqs.(A14) and 
(A15) below). It is then clear that the two functions g^{D) and g<s>(D) cannot be space- independent, as they should, unless, 
either one restricts to the Einstein-de Sitter case, in which case g^(D) = g<s>(D) = 1, or one imposes suitable restrictions on 
the functional form of the initial potential </?i(x). Specifically, it is easy to show that the separability ansatz (A5), (A6) for 
the gravitational and velocity potentials requires the validity of the constraint 

n2 



F 2 (x)=A~ 



Vvsi(x) 



(A9) 



with K a suitable constant. It is interesting to note that a sufficient condition for the validity of this constraint is that 
the second-order invariant of the initial deformation tensor vanishes, /i 2 (ui) = (this is e.g. the case of one-dimensional 
perturbations, which are described by the Zel'dovich solution). This can be easily shown by replacing the relation (A9) into 
the last of Eqs.(26). 



APPENDIX A3: SECOND-ORDER EULERIAN PERTURBATION THEORY IN AN EINSTEIN-DE 
SITTER MODEL 

In this appendix, we want to apply our formalism to the case of an Einstein-de Sitter universe, without cosmological constant. 
In particular, and as an example of how our formalism works, wc want to derive the expressions of <^ 2 '(x, o) and <I>' 2 ^(x, a), 
then the second-order corrections <r 2 ' (x, a) = a 2 82 (x) and u' 2 ' (x, a) = a u 2 (x) already known in the literature. Again, we 
use in this appendix the same notation of the main text. 

The fundamental equations for the cosmological potentials ip and $ are given by 



! + 5<™> a + l;<* + *> 







a 



(aifi) + $ + F = , 



(A10) 



V 2 F = aV • [VVV$] 



As in Section 3, to find approximate perturbative solutions of these equations, we expand ip, $ and F according to the 
relations 



ip{x, a) = <p w (x) + ip (2) (x, a) + 



ifii (x) +a<^ 2 (x) + ■ 



hWl 



s( 2 ) 



$(x,a) = $ ilJ (x) + $ w (x,a) + --- = $!(x)+a$ 2 (x) + - 



(AH) 
(A12) 



and -F(x, a) — F^ 2 '(x, a) = aF2(x); we remind that ^'(x, a) — aVV (x), while <5i(x) = V^i(x), where (T^x, o) = a5i(x); 
<pi(x) = — $i(x). Remind that F is by definition at least a second-order quantity. Finally, we stress that the general expansions 
(24) and (25) have become separable in space and time in the flat case. In terms of the perturbed quantities one has 



16 P. Catelan, F. Lucchin, S. Matarrese and L. Moscardini 
|$ 2 + I(V* 1 ) 2 + | ¥>2 = 0, 

2<^2 + $2 + F 2 = , 

V 2 F 2 = -S? + (ui • V) <Ji . 
From the first two equations it is possible to obtain the expressions for ip 2 and $2: 

1P2 = --F 2 + -(X7<pi) 2 , 

$2 = ^F 2 -^(V^i) 2 . 

We stress again that the velocity field u is manifestly irrotational in physical space, namely 



u(x, a) = V 



?( 2 )( 



$ w (x) + - F^> (x, o) - ^ a u w (x)' 



The final expressions for 8 and V • u are 



<5 (1) +V 2 



aF 



(2) 



1 



<5 (1) (x,a) + ^ (1) (x,a) 2 -a(u (1) • V)<5 (1) (x,a) + 1 a 2 ^(arf) , 



a a 



and 



aV-u = 5 llj +V 



3 r-,(2) 2 2 (1)2 

aF y ' H — a ir ; 

7 7 



= J< 



1} (x, a) + I 8 W (x, a) 2 - a (u (1) • V) S (1) (x, a) + ± a 2 J^ ( d P w « 



a/3 



(A13) 



(A14) 
(A15) 

(A16) 



(A17) 



(A18) 



which give equations (42) and (43) of the main text; they can be also compared with equation (18.8) in Peebles (1980), 
where, however, a slightly different definition of gravitational potential is used, with equations (19a,b) in Wise (1988) and 
with equations (8) and (10) in Catelan & Moscardini (1994b). 



APPENDIX A4: SECOND-ORDER EULERIAN PERTURBATION THEORY IN A CLOSED 
FRIEDMANN MODEL 



The Q > 1 case may be treated starting from equation (44) and transforming to the new variable (see e.g. Peebles 1980) 

(A19) 



-B- 



where now < r < 1, the lower limit corresponding to the Big Bang, the upper limit to the instant of maximum expansion: 
in this appendix we consider only the expanding phase. 

The equation to be integrated is formally identical to equation (46) and the initial conditions are the same. Finally, all 
the non-trivial expressions valid for a closed model are 

D{t) 



J(t) 



D'(t) 



rf(ry 



and J(t) — > t in the limit r — > 0, like in the open models; furthermore 
D(t) 



-l + h^Mir) 



(A20) 



(A21) 



with M(t) = arctan-y/r/(l — r) ; D(t) — 2r/5 in the limit r — > 0, as in the case fi < 1. The function A(t) is now given by 



Mr) = 2 + - 



(l-T)f(r) 2 

3y / r(l-r) - (3 - 2t)M(t) 

(3-TVr(l-T)-(3-3r)M(r) 

with /(0) = 1. We then have: 



and 



/(r) = - 5 



(A22) 



(A23) 
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B(r) 



4D(t 



+ 



4tD(i 



1 — T , ,, s 1 

— M(r)-- 



1 - r M(r) 



(A24) 



If, alternatively, one prefers to work directly in terms of the time-variable f2, the formulae which have to be used are 
respectively: 

D(fl) 



j(n) = 



D'(fi) 



D(Q) = 



m) 



= n(n-i)/(n)- 1 , 

[(1 + 2n)V^ - 1 - 3ft M(Q,)] 



(n-i) a / 2 

Sy 7 ^ - ! - (2 + fi) M(ft) 



3ft 
~2~ 



(1 + 2ft) v 7 ^ - T - 3ft M(ft) 



(A25) 
(A26) 

(A27) 



where M(ft) = arctan^ft - 1 and M(l) = 0. 
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